class Transdim_MCMC
{
	public:
		Transdim_MCMC(mat v,mat X,double *Y, int ds)
		{
			_v=v;
			_ds=ds;
			_X=X;
			_p=X.n_cols-1;
			_n=X.n_rows;
			_Y=Y;
			U=new Distribution::Uniform(1,0,1);
			boost::random::uniform_01<> U;
			runif=new RandomG::Random<boost::random::uniform_01<> >(U);
			fL=new Laplace;
			_d=new boost::math::normal_distribution<>(0,1);
			_is=inv(v);
		}
		vector<mat> Move_theta(mat gamma,mat thetak,mat gammat)
		{
			int nkp=as_scalar(sum(gamma,0));
			int nb=nkp+1;
			vector<mat> vect;
			vector<mat> res;
			mat gammati(_p+1,1);
			mat X2(_n,nb);
//			cout << "test"<< p+1;
			int kk=1;
			for(int i=0;i<(_p+1);i++)
			{
					
				if(i!=0){
					gammati(i,0)=gamma(i-1,0);
					if(gamma(i-1,0)!=0)
					{
						X2(span::all,kk)=_X(span::all,i);
						kk++;
					}

				}else{
					X2(span::all,0)=_X(span::all,0);
					gammati(i,0)=1;
				}
			}

			mat Xt=_X*diagmat(gammati);
			mat v1(nb,nb);
			v1.fill(10);
			mat m(nb,1);
			m.zeros();
			mat v=diagmat(v1);
			(*fL)(X2,_Y,m,_v(0,0),m);			
			mat Skp=(chol((*fL).Get_Sig())).t();	
			mat Mkp=(*fL).Get_Mu();
			mat thetakp;
			double g;
			if(nkp<_nk)
			{
				mat foo=inv(_Bk)*(thetak-_mk);
				thetakp=Mkp+Skp*foo(span(0,nkp),0);
				mat bar=inv(Skp)*(thetakp-Mkp);
				mat u=bar(span(0,_nk-nkp-1),0);	
				mat m(_nk-nkp,1);
				m.fill(0);
				mat s(_nk-nkp,_nk-nkp);
				s.eye();
				Distribution::Gaussian G(_nk-nkp,m,s);
				g=(double)1/G.d(u);

			}else{
				mat m(nkp-_nk,1);
				m.fill(0);
				mat s(nkp-_nk,nkp-_nk);
				s.eye();
				Distribution::Gaussian G(nkp-_nk,m,s);
				mat u=G.r(1);
				mat foo=inv(_Bk)*(thetak-_mk);
				mat bar=rbind(foo,u);
				thetakp=Mkp+Skp*bar;
				g=G.d(u);
			}
			mat beta=add_mat(thetakp,gammati,_p+1);
			double post=exp(Lik_Prob(beta,Xt)+Priorb(thetakp)+Prior(gammati));
			double rho=(det(Skp)*post*g)/(det(_Bk)*_post);
			double un=(*runif)();
			if(un<rho)
			{
				_Bk=Skp;
				_mk=Mkp;
				_nk=nkp;
				_post=post;
				res.push_back(thetakp);
				res.push_back(gamma);
				
			}else{
				res.push_back(thetak);
				res.push_back(gammat);
			}
			return res;
		}
		
		inline double Lik_Prob(mat theta1,mat Xt)
		{
			double L=0;
			double sum=0;
			int m=theta1.n_cols;
			//Xb
			int n=_n;
			double* Y=_Y;
			mat foo(1,Xt.n_cols);
			for(int i=0;i<n;i++)
			{	
				foo= Xt.row(i);
				sum=dot(foo,theta1);
				
				if(Y[i]==1){
					L+=log(Phi(sum));
				}else{
					L+=log(1-Phi(sum));
				}	
			}
			//cout << exp(L);
			return L;

		}
		inline double Phi(double x)
		{
			double t=cdf((*_d),x);
			if(t==1)
			{
				t=0.99999;
			}
			if(t==0)
			{
				t=0.00001;
			}
			return t;
		}
		inline double Priorb(mat theta)
		{
			double res=0;
			mat foo=theta.t();
			//cout << foo << "\n";
			//cout << "//" << bar;
			int d=theta.n_elem;
			mat ss(d,d);
			ss.fill(_v(0,0));
			ss=diagmat(ss);
			double  bar=as_scalar(foo*inv(ss)*foo.t());
			res=-0.5*d*log(2*PI)-0.5*log(det(ss))-0.5*bar;
			return res;

		}
		double Prior(mat theta){
		    
			double t=0;
			double m=0.5;
			double ds=logfact(_ds);
			double d=logfact(_p);
			double i=sum(sum(theta));
			if(i<=_ds)
			{
				double t=logfact(i)+(logfact(_ds-i))-ds+i*log(m)+(_ds-i)*log(1-m)+logfact(i)+(logfact(_p-i))-d;
				return t;
 			}else{
				return -numeric_limits<double>::infinity();

			}		
			
	//		return -0.693;
		}

		
		mat Q_gamma(mat gammat)
		{
			int d=sum(sum(gammat));
			int dc=_p-d;
			 
			double un=(*runif)(); 
		 	if((dc==0) | (dc==_p) | un<(double)1/3)
			{	
				Distribution::Uniform U1(1,0,_p-1);
				int u=as_scalar(U1.r(1));
				gammat(u,0)=1-gammat(u,0);//as_scalar((*U).r(1));	
				return gammat;
			}else{
				
				Distribution::Uniform U1(1,0,d-1);
				Distribution::Uniform U2(1,0,dc-1);
				int u1=as_scalar(U1.r(1));
				int u2=as_scalar(U2.r(1));
				int b=1;
				int i=0;
				int k=-1;
				int j=-1;
				int b1=0;
				int b2=0;
				while(b)
				{
					if(gammat(i,0)==1)
					{
						k++;
					}else{
						j++;
					}
					if(k==u1)
					{
						gammat(i,0)=0;
						b1=1;
					}
					if(j==u2)
					{
						gammat(i,0)=1;
						b2=1;
					}
					i++;
					b=1-b1*b2;
					
				//	cout << gammat;
				}
				
			//	cout << "\\\\\\\\" << gammat;
				return gammat;
			}
		}
		vector<mat> Move(mat gammat, mat thetat)
		{
			mat gamma=Q_gamma(gammat);
			vector<mat> x=Move_theta(gamma,thetat,gammat);
			return x;
		}
		mat Sample(int M)
		{
			mat Gamma(M,_p);
			Gamma.fill(0);
			mat m(_p+1,1);
			m.fill(0);
			Distribution::Binomial U1(1,0.1,_p);
			Gamma.row(0)=U1.r(1).t();
			mat Y=vect2mat<double>(_Y,_n,1);
			double t=as_scalar(mean(Y,0));
			post_init(Gamma.row(0).t());
			mat theta=_mk;
			
			for(int i=1;i<M;i++)
			{
				vector<mat> vect=Move(Gamma.row(i-1).t(),theta);
				Gamma.row(i)=vect[1].t();
				theta=vect[0];
				cout << i << " ";
				cout << vect[1].t();

			}
			return Gamma;

		}
		void post_init(mat gamma)
		{

			int nkp=as_scalar(sum(gamma,0));
			int nb=nkp+1;
			vector<mat> vect;
			vector<mat> res;
			mat gammati(_p+1,1);
			mat X2(_n,nb);
//			cout << "test"<< p+1;
			int kk=1;
			for(int i=0;i<(_p+1);i++)
			{
					
				if(i!=0){
					gammati(i,0)=gamma(i-1,0);
					if(gamma(i-1,0)!=0)
					{
						X2(span::all,kk)=_X(span::all,i);
						kk++;
					}

				}else{
					X2(span::all,0)=_X(span::all,0);
					gammati(i,0)=1;
				}
			}
			cout << gammati;

			mat Xt=_X*diagmat(gammati);
			mat v1(nb,nb);
			mat S(nb,nb);
			v1.fill(10);
			S.fill(1);
			mat m(nb,1);
			m.zeros();
			mat v=diagmat(v1);
			S=diagmat(S);
			(*fL)(X2,_Y,m,_v(0,0),m);			
			_Bk=(chol((*fL).Get_Sig())).t();	
			_mk=(*fL).Get_Mu();
			_nk=nkp;
			G=new Distribution::Gaussian(nb,_mk,S);
			mat res1=((_Bk))*(*G).r(1);
			mat beta=add_mat(_mk,gammati,_p+1);
			cout << beta;
			_post=exp(Lik_Prob(beta,Xt)+Priorb(_mk)+Prior(gammati));

		}

	
	private:
		Distribution::Gaussian *G;
		mat _v;
		mat _is;
		boost::math::normal_distribution<> *_d;
		Laplace *fL;
		Distribution::Uniform *U;
		double *_Y;
		mat _X;
		int _p;
		int _ds;
		int _n;
		mat _Bk;
		int _nk;
		mat _mk;
		double _post;
		RandomG::Random<boost::random::uniform_01<> > *runif;
};

